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Abstract 

We present highly accurate Monte Carlo results for simple cubic Ising lattices 
containing up to 256 3 spins. These results were obtained by means of the Cluster 
Processor, a newly built special-purpose computer for the Wolff cluster simulation 
of the 3D Ising model. We find that the magnetization M(t) is perfectly described 
by M{t) = (ao - ait e - a 2 t)t p , where t = (T c - T)/T c , in a wide temperature range 
0.0005 < t < 0.26. If there exist corrections to scaling with higher powers of t, they 
are very small. The magnetization exponent is determined as (3 = 0.3269 (6). An 
analysis of the magnetization distribution near criticality yields a new determination 
of the critical point: K c = J/ksT c = 0.2216544, with a standard deviation of 3-10~ 7 . 



We consider the 3D Ising model on the simple cubic lattice, with nearest-neighbor 
interactions J, at a temperature T and a coupling strength K = J/kgT. At criti- 
cality, the spontaneous magnetization vanishes with a singularity M(t) oc t 13 , where 
t = (T c — T)/T c = (K — K c )/K parametrizes the distance to the critical point. 
However, this law applies only in the limits of infinite system size and t — > 0. Even 
for the infinite system there are corrections due to nonzero values of t: 

M(t) = f t (t)t? (1) 

where ft is some function of t, finite at t — 0. 

For the 2D Ising model this function is known exactly, and it is analytic. How- 
ever, in the 3D case, f t is not analytic at t — 0. The leading terms of its expansion 
near t = are 

f t (t) « a - a x t e (2) 

where 6 ~ 0.5 is Wegner's correction-to-scaling exponent ]1J. 

Generally one would expect that there exist many more terms, containing higher 
powers of t, in the expansion of ft- Quite remarkably, we find that it is sufficient to 
add only one term, (— c^t), to Eq. (H), in order to describe M{t) of the simple cubic 
3D Ising model with very high accuracy. This does not only apply close to t — 0; it 
holds in a wide range of t. 

We found this intriguing fact using the Cluster Processor (CP) 0. The CP 
implements the cluster Wolff algorithm || in hardware, for 3D simple cubic Ising 
models with nearest-neighbor interactions and periodic boundaries. Its memory and 
speed are sufficient to simulate Ising systems containing up to 256 3 spins. The CP 
was checked to give correct results for 2D Ising systems. Moreover the CP data are 
also consistent with earlier very accurate simulations on 16 3 and 32 3 lattices [Q. In 
the CP, the system size can take the 5 values 16, 32, 64, 128, 256 along each spatial 
direction. But the present work is restricted to systems with the same size L in all 
three directions. 

To determine ft, one should eliminate the influence of the finite system size. This 
can be easily achieved by a comparison of data for different sizes. The corrections 
due to finite L can be described by a scaling function ft,: 

M(t) = f L (Lnt? (3) 

where the exponent v describes the divergence of the correlation length when the 
critical temperature is approached. The relation Eq. (|3]) is valid in the limit t — > 0. 
The argument of the function Jl is proportional to the ratio of the finite size and 
the correlation length. Thus, fi is expected to be a constant for large values of its 
argument. 

The meaning of M(t) in Eq. (|3]) has to be made more precise, since the aver- 
age magnetization of a finite system vanishes. Instead, we may define a nonzero 
expectation value in terms of the absolute value of M: 

M m = (\M\) (4) 

We can also find M(t) from the spin-spin correlation function, using the relation 
M 2 = (S '(0) S '(oo)) . For a finite lattice this infinite distance can be replaced by 
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L/2, half way to the periodic images of the spins. Therefore we use the following 
expression as the second definition of the magnetization: 

M c =(S(0)S(L/2))^ (5) 

Finally, the magnetization can be defined as M 2 = (M 2 ) 1 / 2 . This quantity has been 
studied in for the 2D Ising model and in || for the 3D case. In agreement with 
H || we found that M 2 is strongly affected by finite-lattice effects even for large t, 
so we do not use M 2 in this paper. 

The scaling functions fi are different for M m and M c . This helps to determine 
the ranges of t where the finite-size effects become important. 

Our data show, that for a given L, the finite-size corrections to M m and M c are 
smaller than the error bars for t > t^. For L = 32 this is illustrated by Fig. 1, 
which clearly demonstrates that M m and M c coincide for t > t^2, and £32 is about 
1.5 x 1CT 2 . From Fig. 1 one can see that for t just below ti the corrections to M m 
and M c have different signs, which facilitates determination of ti- 

The value of ti can be estimated in the following way. The finite-lattice ef- 
fects become important when the correlation length £ oc t~~ v is comparable with L. 
Therefore 

at~ L u = L (6) 

where a is a numerical coefficient. Taking into account v = 0.63 Q, we get a « 2, 
which seems reasonable. According to Eq. (Q), a doubling of L decreases by a 
factor 2 1 /" « 3. This is in agreement with Fig. 2, which shows the normalized M m (t) 
results for three different lattice sizes. 

To study f t (t) we used data for lattice sizes ranging from 32 3 to 256 3 . Only 
data for t > t^, where M m and M c coincide, were taken into account. The results 
for M(t) are shown in Fig. 3, using a logarithmic scale on both axes. Close to 
the critical point, for t < 0.02, the plot seems to be linear. However, attempts to 
approximate the data by t /3 P(t), where P(t) is an arbitrary polynomial in t, were 
not successful. Fig. 4 shows the poor result of such an attempt. It describes the 
ratio of the simulation data to 

M int (t) = tP(po + Pl t + p 2 t 2 + p 3 t 3 ) (7) 

The exponent (3 and the coefficients Pi were determined by the least-squares method, 
in order to describe the simulation data as closely as possible by Eq. (|7|). Never- 
theless, the differences between the CP data and the approximation by Eq. (|7|) are 
much larger than the error bars. 

This suggests that we should use the form Eq. (|2|) instead to describe ft(t). Thus 
we wrote f t (t) as a polynomial in t 1 / 2 

fa = o - M 1/2 - a 2 t (8) 

Even the first attempt to approximate the simulation data as t 13 f t i with (3 = 0.3267, 
taken from ||], and a ,ai,a 2 regarded as free parameters, was very successful. But 
the high statistical accuracy of the CP data allows the use of even more adjustable 
parameters, and we supposed that the magnetization can be described by 

M (t) = (a - a x t e - a 2 t)t p (9) 
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in the interval 0.26 > t > 0.0005. The five parameters /?, 9 and a, were determined 
by the least-squares method. To estimate the influence of the uncertainty in K c , we 
fitted the parameters not only for our best estimate K c = 0.2216544 (see below), 
but also for the lower and upper limits of K c , defined by one standard deviation 
(3 x 10~ 7 ) of the critical coupling. The results are shown in Table 1. The last 



Table 1: Results of least-square fits of the magnetization data obtained by the CP. 
These fits were made for three choices of the critical coupling K c . 





P 


9 


a 


a,\ 


a 2 


x 2 


0.2216544 
0.2216541 
0.2216547 


0.3269(3) 
0.3274(3) 
0.3265(3) 


0.508(15) 
0.490(15) 
0.528(15) 


1.692(4) 
1.698(4) 
1.686(4) 


0.344(6) 
0.340(5) 
0.348(7) 


0.426(11) 
0.436(10) 
0.414(12) 


0.844 
0.859 
0.848 



column characterizes the quality of the least-squares approximation. It is defined 
as X 2 — J2iLi ( j^^y^ j > where n = 5 is the number of fitted parameters, 
N = 45 is the number of magnetization data points, and the cr(tj) are the standard 
deviations of the magnetization data M(ti). The minimum of \ 2 as a function 
of K appears to occur near K c = 0.2216544 as determined from an analysis of 
the magnetization distribution near the critical point (see below). This agreement 
between two different approaches suggests that our scaling formulas are adequate. 

^From the data in the second column we estimate the magnetization exponent 
as /3 = 0.3269 (6), which is in a good agreement with earlier values, see e.g. Ref. 0] 
and references therein. 

The third column indicates that 9 = 0.508 (25), supporting earlier results ob- 
tained by means of an e-expansion analysis f7|, series expansions || and a finite-size 
scaling analysis of three different Ising models ||. 

It should be emphasized that the accuracy of Eq. @ is much better than one 
might expect from the quoted standard deviations of the formula parameters. To 
utilize this property, we present here a very accurate empirical approximation for 
the spontaneous magnetization of the simple cubic Ising model in the region 0.26 > 
t > 0.0005: 

M (t) = t a32694109 (1.6919045 - 0.3435773 lt a50842026 - 0.42572366t) (10) 

where t = 1 - 0.2216544 k B T/ J. 

In the range of t between 0.24 and 0.17 the results of Eq. (|T0"D numerically 
coincide with the Pade approximant of Ref. || within 10~ 5 , while for smaller t the 
approximation Eq. ( |i"0"|) is superior. 

The ratio of the CP data to M (t) is shown in Fig. 5. The simulation data 
coincide with M (t) within the error bars. No systematic deviation of the CP data 
from M (t) can be found. 

The function f t may also contain a term proportional to t 26 . Because 9 is very 
close to 0.5, the simulation accuracy is not sufficient to distinguish the term (a^t) 
in Eq. (|]) from the sum (a 2 it + g^ 26 *)- 
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The CP was also used for some additional calculations close to the critical point. 
We obtained the dimensionless ratio Q = (m 2 ) 2 / (m 4 ), related to the Binder cumu- 
lant [ nj . The Q values are shown in Table |2|. 



Table 2: Numerical results for the dimensionless ratio Ql = (^ 2 )|/(^ 4 )l for finite 
three-dimensional Ising models close to the critical point. Also shown is the total 
number of Wolff clusters flipped by the CP to obtain each numerical result. 



L 


K 


Ql 


# clusters 


16 


0.2216530 


0.63381 (5) 


1 x 10 9 


32 


0.2216530 


0.62883 (10) 


3 x 10 8 


64 


0.2216530 


0.62469 (79) 


2.5 x 10 7 


64 


0.2216545 


0.62670 (25) 


2.5 x 10 8 


128 


0.2216530 


0.62189 (168) 


1 x 10 7 


128 


0.2216545 


0.62603 (84) 


4 x 10 7 


256 


0.2216530 


0.61555 (180) 


2.6 x 10 7 



These data were combined with those available from Ref. M in order to determine 
the critical point more accurately. Near T c , the bulk correlation length satisfies 
£ ^> L, so that Lt u , and hence (K — K c )L Vt , is small. Thus Ql{K) can be expanded 

as 

Q L {K) = Q + qi {K- K C )L* + q 2 (K - K C ) 2 L 2 ^ + + b 2 U» (11) 

The renormalization exponents yt and y- x are related to v and 9 as yt=v^ 1 and 
yi=—6/u. The substitutions of the values y t = 1.587(2), y- x = —0.82(6) and y 2 = 
— 1.963(3), taken from ||, into Eq. ([□]) was found to describe the combined data sat- 
isfactorily. A least-squares fit was used to determine K c . For system sizes L > 5 the 
same fit as in Q yielded K c = 0.2216544(3), where the standard error includes the 
uncertainty in the input parameters. As a final estimate we quote K c = 0.2216544(6) 
with an error of two standard deviations, in order to account for a possible bias in- 
troduced by our choice of the form of Eq. (0). The additional data in Table 2 
permitted a clear improvement of the accuracy in comparison with Ref. 

We calculated the irrelevant exponent y- x in two ways. First, it can be obtained 
as a product of 9 = 0.508(25), found from Table (1), and y t =1.587(2) §. The 
result is y\ = —0.81 (4). Another possibility is to include y x as a free parameter 
in the Ql(K) fitting procedure. This yields y\ = —0.83 (9), in agreement with the 
result given above, and with earlier results ||, based on an analysis including three 
different Ising models. Our estimates of y- y agree with e-expansion results of 0. 

Furthermore, we have combined the new data for (m 2 ) near the critical point 
with those obtained in Ref. Q for the nearest-neighbor model with system sizes 
up to L = 40. The analysis was based on the expected scaling behaviour of the 
susceptibility (see Ref. ||) 

L d (m 2 ) = c + Cl (K-K c ) + --- 

+ L 2 y*- d k + di(K - K c )Ly< + d 2 (K - K C )L 2 ^ + 9l L* + g 2 L 2 y' 
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with K c = 0.2216544 and y' = —2.1. The renormalization exponent is related 
to the magnetic susceptibility critical exponent 7: 2y^ — d = 7/1A The analysis for 
L > 5 yielded y^ = 2.4808 (16) where we again quote a two-sigma error. This result 
is in a good agreement with Ref. |4[] and references therein, and with Ref. fTT| . 
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Figure captions 

Fig. 1: Normalized magnetization M m (t) and the correlation function (S(0)S(L/2)) = 
M 2 C for the 32 3 lattice. M (t) is given by Eq.@. The normalization makes it possible 
to expand the scale so that the small deviations of M(t) from Mo(t) become visible. 
(S(0)S(L/2)) is close to M 2 (t), and M(t) is close to M (t) for t > t 32 « 0.015. 
For t > 0.26 the deviations of the CP data from M (t) grow rapidly. In this low- 
temperature range the simulation data are in an excellent agreement with the Pade 
approximant of Ref . || . 

Fig. 2: Normalized magnetization M m (t) for 3 different lattice sizes L = 32, 64, 
and 128. The small-t behavior displays the characteristics of the finite-size scaling 
function fz associated with M m . 

Fig. 3: Spontaneous magnetization of the 3D Ising model as a function of t, with a 
logarithmic scale on both axes. The statistical errors in these data points are indi- 
cated, but they appear only as single horizontal bars, because the errors are below 
the resolution of this figure. The data shown here apply to different lattice sizes, 
and were taken in those ranges of t were finite-lattice effects are negligible: they 
describe the infinite system magnetization. 

Fig. 4: Ratio of magnetization data to the approximation M int (t) given in Eq. ([|), 
where the function f t (t) was supposed to contain only integer powers of t. This Fig- 
ure demonstrates that, without the Wegner correction to scaling, even a 5-parameter 
fit according to Eq. (|7|) does not describe the CP data properly. In this equation we 
used the same number of adjustable parameters as in M (t) (see Eq. |]), but, as can 
be seen by comparing Figs. 4 and 5, the latter approximation is incomparably better. 

Fig. 5: Normalized magnetization data, describing the infinite 3D Ising system. 
This Figure demonstrates that the expression Eq. (H) for M(t) agrees with the sim- 
ulation data within the statistical errors. The relative accuracy of the formula is as 
high as 10 -3 for t ~ 10 -3 , and better than 10 -4 for t > 0.01. This picture combines 
all the magnetization data, obtained with different precisions and for different lat- 
tice sizes. Long simulation times were required for some of these points in order to 
obtain such small error bars. 
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